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The evolution equation for the flame surface 
density in turbulent premixed combustion 

By Arnaud Trouve 


1. Motivation and objectives 

The mean reaction rate in flamelet models for turbulent premixed combustion 
depends on two basic quantities: a mean chemical rate, called the flamelet speed, 
and the flame surface density (Section 1.1). Our previous work had been primarily 
focused on the problem of the structure and topology of turbulent premixed flames, 
and it was then determined that the flamelet speed, when space-averaged, is only 
weakly sensitive to the turbulent flow field (Haworth & Poinsot 1992, Rutland & 
Trouve 1991, Trouve 1991). Consequently, the flame surface density is the key quan- 
tity that conveys most of the effects of the turbulence on the rate of energy release. 
In flamelet models, this quantity is obtained via a modeled transport equation called 
the E-equation. Past theoretical work has produced a rigorous approach that leads 
to an exact but unclosed formulation for the turbulent E-equation (Section 1.2). In 
the exact E-equation, it appears that the dynamical properties of the flame surface 
density are determined by a single parameter, namely the turbulent flame stretch. 
Unfortunately, the turbulent flame stretch as well as the flame surface density is 
not available from experiments, and, in the absence of experimental data, little is 
known on the validity of the closure assumptions used in current flamelet models. 

Direct Numerical Simulation (DNS) is the alternative approach to get basic in- 
formation on these fundamental quantities. In the present work, three-dimensional 
DNS of premixed flames in isotropic turbulent flow is used to estimate the different 
terms appearing in the E-equation (Section 1.3). A new methodology is proposed to 
provide the source and sink terms for the flame surface density, resolved both tem- 
porally and spatially throughout the turbulent flame brush (Section 2.1). Using this 
methodology, our objective is to extract the turbulent flame stretch from the DNS 
data base and then perform extensive comparisons with flamelet models (Section 
2.3). Thanks to the detailed information produced by the DNS-based analysis, it is 
expected that this type of comparison will not only underscore the shortcomings of 
current models, but also suggest ways to improve them. 

1.1 The flamelet approach for turbulent premixed combustion 

Flamelet models constitute one of the most common approaches for turbulent 
premixed combustion (Bray 1980, AVilliams 1985, Peters 1986). Such models are 
very attractive as they conveniently separate the combustion problem from the 
analysis of the turbulent flow field. These models are based on the flame sheet 
assumption, which requires chemical reaction to occur at fast time scales and short 
length scales relative to the turbulence. In this situation, the flame is confined 
to relatively thin layers within the turbulent flow field, and the topology of the 
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reaction zone corresponds to a surface burning mode. Experimental, theoretical, 
and numerical evidence suggest that the domain of validity of the flamelet picture 
is indeed large and that many technologically important reacting flows occur in the 
flamelet regime. 

In the flamelet regime, it is convenient to describe the flame-flow interactions 
in terms of the following two basic ingredients: a flame speed that characterizes 
the flame inner structure and the flame front surface area. For instance, the mean 
reaction rate may be written as a mean integrated chemical rate times the flame 
surface density: 

(ur) = ( j u R dn) s (S'), (1) 

where ljr is the mass of fuel consumed per unit time and per unit volume; f uj R dn 
is the local integral of the reaction rate along the flame normal direction and charac- 
terizes the local combustion intensity; £' is the flame surface area per unit volume. 
The flame surface density is defined as the expected value for £': £ = (£'). Note 
that in Eq.(l), while ( ) denotes a standard, unweighted ensemble-average, ( )$ de- 
notes a surface mean, defined as an area-weighted ensemble- average (see Section 

i- 2) - 

In Eq.(l), the integral of the reaction rate is simply related to the local fuel 
consumption speed, 5c, defined as: 

Sc = — -7 — [ Lj R dn, (2) 

PuYr } uJ 

where p u and Yr }U are respectively the density and the fuel mass fraction in the 
unburnt gas. The mean reaction rate may then be expressed as the surface mean 
of 5c, called the flamelet speed, times the flame surface density: 

(ur) = (PuYr^ u (Sc)s) E. (3) 

Eq.(3) is the classical flamelet expression for the ensemble mean of the reaction 
rate. In general, this mean reaction rate is a function of both position and time, 
(o?/t}(x, t ). An overall reaction rate can also be obtained by space- averaging Eq.(3) 
over an arbitrary volume, V , that encloses the reaction zone: 

Wj(t) = (u R )(x,t)dV = (, p u Y R<u S c (t )) (4) 

where 5\/(£) is the total flame surface area within the volume V, Sy(t) = f £ dV; 
and Sc(t) is a mean fuel consumption speed, 5c(0 = / (Sc) S £dV ' / f EdV. Sc 
can also be expressed as the area-weighted, space- averaged value of Sc integrated 
along the turbulent flame surface: Sc = f Sc & A/ f <5A, where 8 A designates the 
area of an infinitesimal flame element and where integration takes place along the 
flame surface. Sc and Sy are global flame properties and are simply related to the 
turbulent flame speed (see Section 2.2). 
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The flamelet speed, (Sc)s, in Eq.(3) and the mean fuel consumption speed, Sc, 
in Eq.(4) account for local variations of the reaction rate along the flame surface. 
Laminar flame theory indicates that the local flame structure is modified by flow 
divergence, usually characterized by the strain rate acting in the flame tangent 
plane as well as by flame front curvature (Clavin 1985, Williams 1985). Under 
certain conditions, these variations can become critical and lead to partial or total 
quenching of the flame. Recent studies using DNS, however, suggest that quenching 
is a rather unlikely event for turbulent premixed flames (Poinsot et al. 1991). In 
addition, although the local combustion intensity may exhibit large fluctuations 
along the turbulent flame front, particularly for non-unity Lewis number flames, 
and although these fluctuations may in turn result in strong spatial variations of 
the flamelet speed, ( Sc)s , across the turbulent flame brush, DNS suggest that 
these departures from the laminar consumption rate always tend to cancel in the 
mean (Haworth k Poinsot 1992, Rutland k Trouve 1991, Trouve 1991). In the 
simulations, the mean fuel consumption speed, Sc, remain within 10% to 30% from 
the one-dimensional, laminar flame speed value, sl. 

Thus, it appears that in the absence of quenching, the mean fuel consumption 
speed, Sc, is only weakly sensitive to the flow field, and the principle effect of 
turbulence is for the fluctuating velocity field to wrinkle the flame and greatly 
increase its surface area. This phenomenon accounts for most of the increase in the 
overall burning rate due to the turbulence and is represented in Eq.(3) through the 
flame surface density, E. 

1.2 The evolution equation for the flame surface density, E 

In current flamelet models, the flame surface density, E, is obtained via a modeled 
transport equation. This equation was first postulated by Marble k Broadwell 
(1977) based on phenomenological grounds. A more rigorous approach was later 
proposed by Pope (1988) and Candel k Poinsot (1990) who derive an exact balance 
equation for the flame surface- to- volume ratio, E': 

^ + V.XE' = (V.X - nn : VX) E', (5) 

at 

where X is the displacement speed of the flame surface, given by the sum of the fluid 
velocity and the flame propagation speed in the normal direction: X = u + urn; n is 
the unit vector normal to the flame surface; and where we use tensorial notations: 

(nn : VX) = 

The right-hand side of Eq.(5) can also be expressed in terms of flame stretch. 
The flame stretch, Jt, is defined as the rate of change of a Lagrangian flame surface 
element, 6 A: 

k = + X.V(6A). (6) 

dt dt 

A more useful expression for k is in terms of strain rate, flame curvature, and flame 
propagation speed (see for example Candel k Poinsot 1990): 


k = ar + 2wk m , 


( 7 ) 
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where ar is the rate of strain acting in the flame tangent plane: ar = V.u— nn : Vu; 
and k m is the flame surface curvature, as given by the divergence of the flame normal 
direction: 2 k m = V.n. In Eq.(7), positive curvature is chosen convex towards the 
reactants. 

Using Eq.(7), the balance equation for the flame surface- to- volume ratio can be 
re-written as: 

dY* 

+ V.XE' = k £'. (8) 

When ensemble- averaged, this equation yields an exact balance equation for the 
flame surface density (Pope 1988, Cant et al. 1990): 

| + V.(X) s S = (l) s E, (9) 

where the flame surface mean of any quantity Q is given by: {Q)s = (QE') /{S') = 
(QE')/E. Note that surface means are different from standard means. In partic- 
ular, the surface mean of a quantity Q is different from the ensemble mean of Q 
conditioned on being at the flame location (see Section 2.1). 

Eq.(9) can be cast in various forms. For modeling purposes, it is useful to split the 
velocity vector into a mean component and a turbulent fluctuation: u = (U) + u'. 
We can then re-write Eq.(9) as follows: 

— + V.(U)E + V.(u') 5 £ 4- V.(um) 5 £ = ( ar)s £ + (At)s £ 4- 2 (ivk m )s £, (10) 
where we use the following notations: 

{ar)s = (V.u ; - nn : Vu') 5 , 

(At}s = ^(U) — (nn ) 5 : V(U). 

The three convective terms on the left-hand side of Eq.(10) are transport terms 
that correspond respectively to convection by the mean flow, turbulent diffusion, 
and flame propagation. The terms on the right-hand side of the equation are the 
source and sink terms for the flame surface density: (clt)s is the turbulent strain 
rate acting in the flame tangent plane, { At)s is the strain rate due to the mean 
flow field, and 2 (wk m }s is a term that accounts for the combined effects of flame 
curvature and flame propagation. 

Eq.(10) is an exact expression of the evolution equation for the flame surface 
density, £, in a turbulent flow field. Unfortunately, many terms remain unclosed in 
this exact turbulent £-equation. Modeling assumptions are, therefore, required, in 
particular to calculate the turbulent diffusion term and the turbulent flame stretch. 
Different formulations of the modeled E-equation have been proposed in the litera- 
ture (Marble &; Broadwell 1977, Darabiha et al. 1987, Maistret et al. 1989, Candel 
et al. 1990, Cant et al. 1990, Borghi 1990, Mantel & Borghi 1991). However, since 
the turbulent flame stretch as well as the flame surface density is not available 
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from experiments, little is known on the validity of the closure assumptions used in 
current flamelet models. DNS is an attractive alternate approach to get basic in- 
formation on these fundamental quantities. In the present work, three-dimensional 
DNS of premixed flames in isotropic turbulent flow is used to estimate the differ- 
ent terms appearing in the E-equation, with particular emphasis on the source and 
sink terms on its right-hand side that determines production or dissipation of flame 
surface area. 

1.3 Direct Numerical Simulation of turbulent premixed flames 
We use DNS to analyze the different terms in the equation for the flame surface 
density. The simulations are performed using a three-dimensional, compressible 
Navier-Stokes solver that fully resolves the turbulent flow field. Spatial derivatives 
are computed with a modified Pade scheme that is sixth-order accurate (Lele 1992). 
Solutions are advanced in time using a third-order Runge-Kutta method (Wray 
1990). Boundary conditions are specified using the NSCBC method (Poinsot & 
Lele 1992). Because of the otherwise prohibitive computational cost, simulations 
are limited to simple but finite-rate reaction schemes. In this work, the chemistry 
model is a single step, irreversible chemical reaction where the reaction rate depends 
exponentially on temperature (Arrhenius kinetics): 

u R - BpY R exp ’ ( n ) 

where T a is the activation temperature and B a constant that depends on the 
selected flame speed. This formulation corresponds to a binary reaction in which one 
of the reactants, Yr , is strongly deficient as, for example, in fuel-lean combustion. 
Also, it is worth emphasizing that the simulations are not limited by the constant 
density assumption, and heat release effects are fully accounted for. 

Following Williams (1985), we re-write the reaction rate as: 

w„=A^e,p( 1 -_^ 1 -.% ), (12) 

where 0 is the reduced temperature, 0 = (T— T u )/(Tt — T„); T u is the temperature 
of the fresh reactants; T b is the adiabatic flame temperature; and the coefficients 
A, a, and /? are, respectively, the reduced pre-exponential factor, the heat release 
factor, and the reduced activation energy: 

A = £exp(— /?/<>), a = (T k -T u )/T t , and fl = aT a /T b . (13) 

In the following, we use a = 0.75 and {3 = 8. 

Another important feature of the simulations is that the transport coefficients 
are temperature dependent. These coefficients satisfy the following relations: 

p = p u (T/T u ) h , Le = \/pDc p = constant, Pr = pc p /\ = constant, (14) 
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where fi, A, and D are the molecular diffusivities of, respectively, momentum, inter- 
nal energy, and species mass, 6 is a constant, and Le and Pr axe, respectively, the 
Lewis number and the Prandtl number. We use b — 0.76, Pr — 0.75. Simulations 
have been performed for different Lewis numbers, Le = 0.8, 1.0, and 1.2. 

The selected computational configuration corresponds to a premixed flame em- 
bedded in a three-dimensional, decaying, isotropic turbulent flow. The left- and 
right-hand sides of the computational domain are inflow and outflow boundaries, 
while periodic boundary conditions are applied at lateral walls. The calculations 
are initialized with fresh reactants on one side of the domain and burnt products on 
the other side; the two are separated by a plane laminar flame. Isotropic turbulence 
is initially located in the flow of fresh reactants, its velocity field being specified 
according to a model spectrum. The turbulence is characterized by a Kolmogorov 
length scale smaller than the thermal thickness of the laminar flame, tj^/St = 0.1 
where St = (Tb — T u )/(dT/dx ) max , and a turbulence intensity that is much higher 
than the laminar flame speed, u f / si — 10. The initial turbulent Reynolds number, 
based on the Taylor microscale, is 50. The initial turbulent Reynolds number, based 
on the integral length scale, is 70. The grid resolution is 129 3 . 

The simulations describe the wrinkling of the flame zone due to turbulent motions 
as well as the combustion feedback due to dilatation and temperature-dependent 
transport properties. Note that the turbulence is decaying in time, and flow condi- 
tions are non-stationary. 


2. Accomplishments 

2.1 Diagnostics 

All terms appearing in Eq.(10) can be obtained from the simulations. We now 
briefly describe how. The velocity vector and the velocity gradient tensor are readily 
obtained from the resolved flow field. To define flame-based quantities, we make 
use of concepts based on the thin flame picture. First, a progress variable, c, is 
introduced to indicate location within the reaction zone, c = 1 — Y/j, where Yr 
is the normalized fuel mass fraction. The progress variable varies monotonically 
through the flame from 0 in the reactants to 1 in the products. Constant progress 
variable surfaces may conveniently be used to define the flame front location: we 
use the surface c = Cf = 0.8. In addition, at any location on this surface, the local 
gradient of c defines the normal direction to the flame front: 


n = — 


Vc 
|Vc| ’ 


(15) 


where n points into the fresh reactants. 

The propagation speed of the flame surface, w , is obtained from an expression 
analog to the well-known field equation (also called the G-equation, see for example 
Kerstein et al. 1988). Let us first consider a point on the flame surface, c — cj. The 
velocity, X, at which this point must move to remain on the surface, is given by: 


Oc 

dt 


+ X.Vc = 0, 


(16) 
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which, using Eq.(15), implies that: 


X.n = 


1 dc 
| Vc| dt ’ 


(17) 


and yields the following expression for the flame propagation speed: 


w 


— X.n - u.n — 


1 Dc 

IVcIDt* 


(18) 


where the quantities axe estimated at the surface c = c/. In the simulations, we solve 
for a conservation equation for the fuel mass fraction, Y/t, and since by definition 
c = x — y R4) the expression above can be readily obtained from the DNS data base. 

The flame surface density, E, is a more subtle quantity. It includes both geomet- 
rical and statistical information. Following Pope (1990), the flame surface density 
is expressed as the product of the expected value for the magnitude of the gradient 
of c, conditioned on being on the flame surface, times the probability of being on 
that surface: 


S' = |Vc| S(c - c/), and E = (£') = <|Vc| 1 c = c f ) p(c/), (19) 

where p(c f ) is the probability of c = c/. An equivalent expression for S, also 
proposed by Pope (1990), is in terms of the joint probability of the scalar c and the 
magnitude of its gradient: 

'' J " /» OO 

S = / 7P( c /i 7)^7> ( 20 ) 

Jo 

where 7 = |Vc|. 

We now turn to the averaging problem. In the simulations, the flame brush 
propagates along the x direction, and the problem remains homogeneous in the 
y — z planes. Therefore, averaged quantities depend on x and time t only, and 
ensemble- averaging can be performed in the y — z planes: 

(Q)(x,t) = j Q{x,y,z,t)dydz, (21) 

where L y and L z are the y and 2 dimensions of the computational domain. Since the 
flow field is non-stationary, time-averaging cannot be used and the present analysis 
is based exclusively on the slicing of the DNS data base in planes perpendicular to 
the direction of mean flame propagation. 

Let us now consider a particular y — z plane, x = xq. In that plane, the flame 
contour c = c/ is a set of closed lines. The flame surface density as well as the 
various surface and conditional means in the analysis can be expressed in terms of 
line-integrals along this (c = Cf,x = xq ) line-contour. For instance, as suggested 
by Eq.(20), the flame surface density, E, can be obtained via the joint probability 
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FIGURE 1. A geometrical expression for the joint probability of the progress 
variable and the magnitude of its gradient in a given y — z plane: p(c/, 7 ) dc d'y = 
prob(cj < c < Cf + dc\ 7 < |Vc| < 7 + cfy) = SS / L y L z , see Eq.( 22 ). 


p(c/, 7 ). Simple geometrical considerations lead to the following relation (Figure 

1 ): 

, J 


p( c f,j)dj = 


E 


dl 


LyLz “ 7 l cos6\' 


( 22 ) 


where dl designates the arc-length of an infinitesimal flame element on the line- 
contour c = Cf and 0 is a measure of the orientation of this flame element with 
respect to the x = plane. In Eq.(22), summation is over all occurrences of the 
event (c, |Vc|) = (c/, 7 ). Eqs.(20) and ( 22 ) then give the following expression for £: 


dl _ 1 

C/ \cos6\ L y L z 


where we use the geometrical relation: |cos#| = y (f ^) 2 + (ff ) 2 /l^ c l- 

Eq.(23) differs from the original formulation proposed by Trouve & Poinsot (1992) 
due to an error in their original expressions for surface means and conditional means. 
However, their estimate of the probability p(c/) remains unchanged: 




P(C/) i P(C/ ’ 7)<f7 L y L z f c=Ct \cos0\\Vc\ L y L z J c=Cf y ( |£ )2 + ( | f )r 

I(24) 

and, compared to the new formulation, the results of Trouve & Poinsot (1992) still 
hold within a few percents. 

Using Eq.( 22 ), the correct expression for conditional means is as follows: 


(Q I c = Cf)(x,t) = 


/ C=C/ Q (|cos0| |Vc |) _1 dl 
J c=c/ (|co^||Vc|)-i dl ■ 


(25) 
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, ,, _ <<?£') _ (QIVcl 1 c_ = cj) _ Wgjf^L 

{Q)s{x,t) (E>) <| Vc || c = C/ ) J c=c/ \cos6\-'dl ' 


(26) 


Clearly, surface means differ from conditional means. 

The relations above provide a methodology to estimate the different means needed 
in our analysis. The accuracy of these expressions depends on the size of the com- 
putational domain with respect to the turbulent length scales. In the simulations, 
the integral length scale of the turbulent flow field grows as the kinetic energy de- 
cays; this growth, however, is rather slow, and it has been determined that the 
integral length scale remains at least 8 times smaller than L y and L x . Although 
we recognize that the statistics axe somewhat undersampled, reasonable accuracy 
is expected when estimating the first moments. 

2.2 Global flame properties 

As described in the previous section, the analysis takes advantage of the fact that 
the problem is statistically one-dimensional and provides the various quantities of 
interest as a function of time t and position x within the turbulent flame brush. We 
now discuss briefly the relationship between these spatially-resolved quantities and 
the more familiar, global (space- averaged) flame properties. 

By definition, the total flame surface area, Sy, within the computational domain 
is equal to the volume-integral of the flame surface density, which, in the present 
case, translates to the following relation: 

= / S (x,t)dx. (27) 

Using Eq.(9), it can be shown that the time-derivative of the relative increase in 
flame surface area is equal to the total flame stretch, k, space-averaged across the 
turbulent flame brush: 

k(t) = = j <*)s(*,0 SOMMr. (28) 


k gives the instantaneous rate of change of the total flame surface area in the 
computational domain. If positive, the flame surface grows; if negative, the flame 
surface contracts. 

The mean fuel consumption speed, Sc > can also be expressed in terms of space- 
integrals of the flamelet speed, ( Sc)s , and the flame surface density, S, across the 
turbulent flame brush: 


Sc(t) = 


f (JtcjsQcjQ s(x,t)di 

f E(x, t)dx 


(29) 
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Finally, using Eq.(4), the turbulent flame speed, St, is simply related to the 
product of Sy and Sc* 

S T (t) = f (u R )(x,t)dx = f (S c ) six, t) S (x,t)dx = Sc(t) (30) 

PuYr.uJ J L y L z 

As described in the previous sections, the present study uses three-dimensional 
direct numerical simulations of turbulent premixed flames in isotropic flow. Three 
different cases have been studied that correspond to turbulent flames characterized 
by the same laminar thermal thickness, <$x> the same laminar flame speed, 
embedded in the same initial turbulent flow field, but with different Lewis number, 
Le = 0.8, 1.0, and 1.2. 

The three cases exhibit large differences in the time history of the turbulent flame 
speed, 5Y(0- As suggested by Eq.(30), these differences may be related to a modi- 
fication of the turbulent flame structure through changes in Sc, to a modification 
of the total flame surface area through changes in Sv , or to a combination of both. 
In fact, irrespective of the value of the Lewis number and in agreement with the 
findings from previous studies (Ashurst et al. 1987, Haworth & Poinsot 1992, Rut- 
land & Trouve 1991), the mean fuel consumption speed, Sc, is found to be only a 
weak function of the turbulence. At all times, departures of Sc from the laminar 
value, $£,, remain within 20%. This result implies that the large effect of the Lewis 
number on St observed in the simulations is primarily a flame surface effect and 
that the flame surface density is the key quantity that to first order determines the 
mean reaction rate. For more details, see Trouve & Poinsot (1992). 

2.S Comparison of DNS results with flamelet models 

Using the methodology described in Section 2.1, we now analyze the spatial distri- 
bution of flame surface density across the turbulent flame brush. Figure 2 compares 
several E-profiles taken at different instants in the simulations. The Lewis number is 
0.8. At t — 0, E is a delta function located at x = 0. As time evolves, the turbulent 
flame brush gets thicker and propagates deeper into the reactants. Accordingly, the 
E-profile spreads out and shifts towards negative values of x. In the simulations, 
this shift is rather weak but can clearly be seen at the latest times ( t = 6 r in Figure 

2 )- 

Figure 3 compares the E-profiles obtained for different Lewis numbers. The 
comparison is performed at the same instant in the simulations, t = 4r. The effect 
of the Lewis number is clearly seen in this plot with a turbulent flame brush that 
propagates faster and deeper into the reactants for more diffusive mixtures. Note 
that as shown in Eq.(27), the z-integral of E gives the relative increase of total 
flame surface area and the local differences in the profiles plotted in Figure 3, when 
integrated, correspond to large global differences in the total amount of flame surface 
area, Sy. 

A different perspective is adopted in Figure 4 where x-location within the tur- 
bulent flame brush is indicated by the local mean progress variable, (c). Although 
they are different in physical space, the E-profiles obtained for different Lewis num- 
bers collapse when mapped in (c) -space. This result indicates that the effect of the 
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x location 

FIGURE 2. Time evolution of the flame surface density, E, in the computational 
domain (reactants on the left; products on the right). Le = 0.8. E and x are made 
non-dimensional by the laminar thermal thickness, St- Time is measured in units 
of the initial turbulent eddy turn-over time, r. 



x location 

FIGURE 3. Lewis number effects on the spatial distribution of flame surface 
density, E, in the computational domain (reactants on the left; products on the 
right), t = 4 t . E and x are made non-dimensional by the laminar thermal thickness, 
St- 

Lewis number on the spatial distribution of flame surface density is primarily an 
effect on the turbulent flame thickness. 

Also plotted in Figure 4 is the flame surface density predicted by the Bray- 
Moss-Libby (BML) model. The BML model is a well-documented, second order 
closure model designed to avoid gradient transport assumptions (see for example 
Bray 1980, and 1990, Libby k Bray 1980). The model pre-assumes a bimodal shape 
for the probability density function of the progress variable, an assumption that is 
consistent with a laminar flamelet description of the reaction zone and represents the 
intermittent presence of packets of unburnt and fully burnt mixture at any location 
within the turbulent flame brush. Various degrees of refinement are possible within 
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FIGURE 4. The spatial distribution of flame surface density, E, across the turbulent 
flame brush: comparison between DNS and the Bray-Moss-Libby model, t = 4r. E 
is made non-dimensional by the laminar thermal thickness, 8 ? . 


the BML framework, in particular with respect to the chemical source terms and the 
flame surface density. Although a transport equation for E could be incorporated 
easily (Cant et aL 1990), current formulations of the model are limited to a simple 
algebraic closure (Bray et aL 1984, 1989). The expression for E in a y — z plane, 
x = xq, where (c) = constant, is: 


E = 


g(c)(l - (c)) 

a yLy 


( 31 ) 


where g is a constant of order unity, L y is a length that characterizes the flame 
wrinkling, and ~a y is the mean of a direction cosine defining the flamelet orientation 
relative to the plane, x = x 0 . Both L y and a y are assumed to be constant across 
the turbulent flame brush. L y is specified empirically as follows: 

Ly = C L l (32) 

where l is the integral length scale of the turbulence, and Cl and n are model 
constants. Following Bray (1990), we use n = 1 , Cl = 1 , g = 1.5, a y = 0.5. 

Figure 4 shows that the structure of the E-profile is not well captured by the BML 
model which predicts a symmetric distribution of E in (c)- space, whereas DNS gives 
a skewed profile with a maximum close to (c) = 0.7. Furthermore, although the 
comparison between DNS and the BML model is rendered difficult by the turbu- 
lence decay, some qualitative observations can still be made about the accuracy or 
inaccuracy of the model. For instance, the BML predictions always overestimate 
the flame surface densities across the turbulent flame brush by a factor decreasing 
from more than 5 to approximately 2 in the course of the simulations. Although 
the exact magnitude of the error might not be meaningful, the fact that it is time- 
dependent underscores the limits of an algebraic closure that cannot account for 
memory effects and is, therefore, expected to perform poorly in transient problems. 
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Figure 5. The spatial distribution of turbulent flame stretch, { k)s , across the 
turbulent flame brush* comparison between DNS and the Coherent Flame Model. 
t — 4 r. k is made non-dimensional by the laminar flame time, 6x/si. 

We now turn to the more promising approach based on a modeled turbulent 
E-equation. Different formulations of this equation have been proposed in the liter- 
ature (Marble & Broadwell 1977, Darabiha et al. 1987, ^4aistret et al. 1989, Candel 
et al. 1990, Cant et al. 1990, Borghi 1990, Mantel &: Borghi 1991). These models are 
based on the exact but unclosed formulation given in Eq.(10). The right-hand side 
of Eq.(10) is the total flame stretch, which can be decomposed into a strain term 
and a propagation term. While the strain term, ( ax)s , remains approximately 
constant through the turbulent flame, the propagation term, 2 {wk m )s, exhibits 
strong variations and decreases from positive values on the unburnt side, close to 
(c) = 0, to negative values on the burnt side, close to (c) = 1 (Trouve & Poinsot 
1992). Interestingly enough, both contributions to stretch have the same order of 
magnitude. 

Figure 5 presents the mean flame stretch profile, ( k)s , across the turbulent re- 
action zone, plotted for different Lewis numbers at t = 4 r. At the leading edge, 
strain rate and flame propagation effects are locally cumulative, and the overall 
balance is strongly positive. The leading edge of the turbulent flame is a region 
of strong production of flame surface area. On the contrary, the propagation term 
takes large negative values on the burnt side. In that region, strain rate and flame 
propagation effects are locally opposite, and the overall balance is negative. The 
rear edge of the turbulent flame thus appears as a region where flame surface area 
gets strongly dissipated. In addition, while the effect of the Lewis number is not 
visible on the burnt side of the turbulent flame, there are large differences between 
the cases Le = 0.8, 1.0, and 1.2 close to (c) = 0. This suggests that the turbulent 
flame is most sensitive to Lewis number effects on the unburnt side, whereas it 
remains relatively unaffected on the burnt side. 

Using Eq. (28), one can quantify the overall differences due to the Lewis number; 
we find (in units of the laminar flame time): k = 2.6 for Le = 0.8, and the flame 
surface area is quickly growing; k — 0.0 for Le = 1.0 and Le — 1.2, and the flame 
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surface area remains roughly constant. 

Also plotted in Figure 5 is a comparison with the Coherent Flame Model proposed 
originally by Marble Broadwell (1977) and further developed by Candel and his 
co-workers (Darabiha et al. 1987, Maistret et al. 1989, Candel et al. 1990). In this 
formulation, the turbulent flame stretch is written as: 

{k)s = vr SL wv (33) 

where k t is the turbulent kinetic energy and e its dissipation, both measured in the 
unburnt gas; {Y r) is the ensemble- averaged fuel mass fraction. The first term on 
the right-hand side of Eq.(33) represents the straining due to the flow motions and 
is assumed to scale with the integral time scale of the turbulence; the second term is 
a disparition term associated with flame propagation and is assumed to scale with 
the laminar flame speed, sl, and the flame surface density, E. 

Figure 5 shows that this model is indeed able to reproduce qualitatively the 
spatial structure of the balance between production and dissipation of E, going from 
production at the leading edge of the turbulent flame to dissipation at the rear edge. 
However, the 1/ (Yr) behavior of the disparition term leads to numerical difficulties 
on the burnt side of the flame. The model, therefore, overpredicts the dissipation 
of E near (c) = 1 and gives a negative mean flame stretch, k = —3.6, in strong 
disagreement with the values reported above. In addition, the disparition term 
in Eq.(33) is always and everywhere negative and cannot account for the possible 
transition to unstable thermo-diffusive conditions as observed in the simulations 
(Trouve & Poinsot 1992). 

3. Conclusion and future plans 

In flamelet models, the flame surface density is a central ingredient that conveys 
most of the effects of the turbulence on the rate of energy release. The flame 
surface density can be obtained using an algebraic closure, as in current formulations 
of the Bray-Moss- Libby model, or via a modeled transport equation, called the 
E-equation. The E-equation was first postulated by Marble & Broadwell (1977) 
based on phenomenological grounds. Recent theoretical work has produced a more 
rigorous approach that leads to an exact but unclosed formulation for the turbulent 
E-equation (Pope 1988, Candel &; Poinsot 1990). Since the flame surface density and 
the different terms of interest in the turbulent E-equation are not experimentally 
accessible, little is known on the validity of the closure assumptions used in current 
formulations of flamelet models. 

Direct Numerical Simulation (DNS) is an attractive alternative approach to get 
basic information on these fundamental quantities. In the present work, three- 
dimensional DNS of premixed flames propagating in isotropic turbulent flow have 
been used to estimate the different terms appearing in the E-equation. A new 
methodology has been proposed to provide the flame surface density, E, as well as 
the source and sink terms for E, estimated as a function of time and position within 
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the turbulent flame brush. Using this methodology, comparisons with flamelet 
models are possible at a level of detail never reached before. 

Principal findings are that: (1) the Bray-Moss-Libby formulation overestimates 
the flame surface densities; more importantly, the error is a function of time, which 
underscores the limits of using an algebraic closure; (2) the Coherent Flame Model 
(CFM), based on a modeled turbulent E-equation, is shown capable of reproducing 
qualitatively the spatial variations in the balance between production and dissi- 
pation of flame surface area — the leading edge of the turbulent flame is correctly 
described as a region of production of flame surface area, whereas the rear edge is 
a region where flame surface gets strongly dissipated - however, the CFM model 
tends to overpredict flame surface dissipation and is unable to reproduce variations 
due to thermo-diffusive phenomena. 

Ongoing plans include the development of a new formulation of the S-equation 
that would incorporate thermo-diffusive mechanisms. Also, the effect of the mean 
flow field as represented by ( At)s in Eq.(10) is currently studied using three- 
dimensional DNS of a turbulent flame embedded in homogeneous shear flow. 

Finally, it is worth emphasizing that theoretical descriptions of turbulent com- 
bustion are aimed at providing tractable expressions for the following two categories 
of unclosed terms: (1) the mean chemical source term; (2) the double and/or triple 
correlations representing turbulent transport phenomena. We have focused exclu- 
sively so far on the first category of problems, i.e. predicting the mean reaction 
rate. The second category is also a promising field for DNS of turbulent combus- 
tion. Indeed, in most current models, turbulent transport terms are treated using 
conventional, first-order closure assumptions based on the concept of an eddy vis- 
cosity. However, since flames of practical interest feature large changes in mass 
density which affect significantly the turbulent flow field, it is somewhat suspicious, 
although very common, to directly extend the empirical closures developed in the 
context of constant density, non-reacting flow problems to the case of turbulent 
flames. Future work will include studying and modeling of turbulent transport in 
a combustion environment. 
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